Go back to the About page.
Let us set some global options for all code chunks in this document.
knitr::opts_chunk$set(
message = FALSE, # Disable messages printed by R code chunks
warning = FALSE, # Disable warnings printed by R code chunks
echo = TRUE, # Show R code within code chunks in output
include = TRUE, # Include both R code and its results in output
eval = TRUE, # Evaluate R code chunks
cache = FALSE, # Enable caching of R code chunks for faster rendering
fig.align = "center",
#out.width = "100%",
retina = 2,
error = TRUE,
collapse = FALSE
)
rm(list = ls())
set.seed(1982)# Install R-INLA package
# install.packages("INLA",repos = c(getOption("repos"),INLA ="https://inla.r-inla-download.org/R/testing"), dep = TRUE)
# Update R-INLA package
# inla.upgrade(testing = TRUE)
# Install inlabru package
# remotes::install_github("inlabru-org/inlabru", ref = "devel")
# Install rSPDE package
# remotes::install_github("davidbolin/rspde", ref = "devel")
# Install MetricGraph package
# remotes::install_github("davidbolin/metricgraph", ref = "devel")
library(INLA)
library(inlabru)
library(rSPDE)
library(MetricGraph)
library(Matrix)
library(dplyr)
library(plotly)
library(scales)
library(patchwork)
library(tidyr)
library(ggplot2)
library(reshape2)
library(sf)
library(here)
library(rmarkdown)
library(knitr)
library(grateful) # Cite all loaded packages
library(latex2exp)
library(plotrix)
rm(list = ls()) # Clear the workspace
set.seed(1982) # Set seed for reproducibilityFollow this link for more details.
# Matern covariance function
matern.covariance <- function(h, kappa, nu, sigma) {
if (nu == 1 / 2) {
C <- sigma^2 * exp(-kappa * abs(h))
} else {
C <- (sigma^2 / (2^(nu - 1) * gamma(nu))) *
((kappa * abs(h))^nu) * besselK(kappa * abs(h), nu)
}
C[h == 0] <- sigma^2
return(as.matrix(C))
}
# Folded.matern.covariance.1d
folded.matern.covariance.1d.local <- function(x, kappa, nu, sigma, L = 1, N = 10, boundary = c("neumann",
"dirichlet", "periodic")) {
boundary <- tolower(boundary[1])
if (!(boundary %in% c("neumann", "dirichlet", "periodic"))) {
stop("The possible boundary conditions are 'neumann',
'dirichlet' or 'periodic'!")
}
addi = t(outer(x, x, "+"))
diff = t(outer(x, x, "-"))
s1 <- sapply(-N:N, function(j) { # s1 is a matrix of size length(h)x(2N+1)
diff + 2 * j * L
})
s2 <- sapply(-N:N, function(j) {
addi + 2 * j * L
})
if (boundary == "neumann") {
C <- rowSums(matern.covariance(h = s1, kappa = kappa,
nu = nu, sigma = sigma) +
matern.covariance(h = s2, kappa = kappa,
nu = nu, sigma = sigma))
} else if (boundary == "dirichlet") {
C <- rowSums(matern.covariance(h = s1, kappa = kappa,
nu = nu, sigma = sigma) -
matern.covariance(h = s2, kappa = kappa,
nu = nu, sigma = sigma))
} else {
C <- rowSums(matern.covariance(h = s1,
kappa = kappa, nu = nu, sigma = sigma))
}
return(matrix(C, nrow = length(x)))
}
# Function to get the true covariance matrix
gets_true_cov_mat = function(graph, kappa, nu, sigma, N, boundary){
h = c(0,graph$get_edge_lengths()[1]*graph$mesh$PtE[,2])
true_cov_mat = folded.matern.covariance.1d.local(x = h, kappa = kappa, nu = nu, sigma = sigma, N = N, boundary = boundary)
return(true_cov_mat)
}
# Define the graph
r <- 1/(pi)
theta <- seq(from = -pi, to = pi, length.out = 100)
edge <- cbind(1+r+r*cos(theta), r*sin(theta))
edges <- list(edge)
graph <- metric_graph$new(edges = edges)
# parameters
h.ok <- 2^-10
type <- "covariance"
type_rational_approximation = "brasil"
rho <- 0.5
#m = 4
sigma <- 1
N.folded <- 10
boundary <- "periodic" # Do not change this
# Mesh sizes
h_aux <- seq(5.5, 4, by = -1/2)
h_vector <- 2^-h_aux
h_label <- paste0("2^-", h_aux, "")
h_label_latex <- sprintf("$2^{-%f}$", h_aux)
# Beta values
beta_aux <- c(4, 6, 6.5, 7, 7.5)
beta_vector <- beta_aux/8
beta_label <- paste0(beta_aux, "/8")
theoretical_rate <- pmin(4*beta_vector-1/2,2)
graph.ok <- graph$clone()
# Build graph with overkill mesh
graph.ok$build_mesh(h = h.ok)
graph.ok$compute_fem()
# Get the overkill mesh locations
loc.ok <- graph.ok$mesh$VtE # or graph.ok$get_mesh_locations()
# Initialize the list of graphs and the list of projection matrices
graphs <- list()
A <- list()
for(i in 1:length(h_vector)){
graphs[[i]] <- graph$clone()
graphs[[i]]$build_mesh(h = h_vector[i])
A[[i]] <- graphs[[i]]$fem_basis(loc.ok)
}
cov.error.folded <- matrix(NA, nrow = length(h_vector), ncol = length(beta_vector))
m_values <- c()
for (j in 1:length(beta_vector)) {
beta <- beta_vector[j]
alpha <- 2*beta
fract2beta <- alpha - floor(alpha)
nu <- alpha - 0.5
kappa <- sqrt(8*nu)/rho
tau <- sqrt(gamma(nu) / (sigma^2 * kappa^(2*nu) * (4*pi)^(1/2) * gamma(nu + 1/2))) #sigma = 1, d = 1
Sigma.folded <- gets_true_cov_mat(graph = graph.ok,
kappa = kappa,
nu = nu,
sigma = sigma,
N = N.folded,
boundary = boundary)
for (i in 1:length(h_vector)) {
h <- h_vector[i]
m <- min(20, ceiling((min(4*beta - 1/2,2) + 1/2)^2*log(h)^2/(4*pi^2*fract2beta)))
m_values <- c(m_values, m)
Sigma <- matern.operators(alpha = alpha,
kappa = kappa,
tau = tau,
m = m,
graph = graphs[[i]],
type = type,
type_rational_approximation = type_rational_approximation)$covariance_mesh()
Sigma.approx <- A[[i]]%*%Sigma%*%t(A[[i]])
cov.error.folded[i,j] <- sqrt(as.double(t(graph.ok$mesh$weights)%*%(Sigma.folded - Sigma.approx)^2%*%graph.ok$mesh$weights))
}
}
print(m_values)## [1] 20 20 20 20 5 4 4 3 4 4 3 2 4 3 3 2 3 3 2 2
slope_circle <- numeric(length(beta_vector))
for (u in 1:length(beta_vector)) {
slope_circle[u] <- coef(lm(log(cov.error.folded[,u]) ~ log(h_vector)))[2]
}
transposed_df <- data.frame(t(data.frame(beta = beta_label, theoretical_rate = theoretical_rate, slope = slope_circle)))
rownames(transposed_df) <- c("beta", "Theoretical rates", "Circle graph")
colnames(transposed_df) <- NULL
# Display the transposed data frame
transposed_df |> paged_table()loglog_line_equation <- function(x1, y1, slope) {
b <- log10(y1 / (x1 ^ slope))
function(x) {
(x ^ slope) * (10 ^ b)
}
}
guiding_lines <- matrix(NA, nrow = length(h_vector), ncol = length(beta_vector))
for (j in 1:length(beta_vector)) {
guiding_lines_aux <- matrix(NA, nrow = length(h_vector), ncol = length(h_vector))
for(k in 1:length(h_vector)){
point_x1 <- h_vector[k]
point_y1 <- cov.error.folded[k, j]
slope <- theoretical_rate[j]
line <- loglog_line_equation(x1 = point_x1, y1 = point_y1, slope = slope)
guiding_lines_aux[,k] <- line(h_vector)
}
guiding_lines[,j] <- apply(guiding_lines_aux, 1, mean)
}
# Generate default ggplot2 colors
default_colors <- scales::hue_pal()(ncol(guiding_lines))
# Create the plot_lines list with different colors for each line
plot_lines <- lapply(1:ncol(guiding_lines), function(i) {
geom_line(data = data.frame(x = h_vector, y = guiding_lines[, i]),
aes(x = x, y = y), color = default_colors[i], linetype = "dashed", show.legend = FALSE)
})
df <- as.data.frame(cbind(h_vector, cov.error.folded))
colnames(df) <- c("h_vector", beta_label)
df_melted <- melt(df, id.vars = "h_vector", variable.name = "column", value.name = "value")
p2 <- ggplot() +
geom_line(data = df_melted, aes(x = h_vector, y = value, color = column)) +
geom_point(data = df_melted, aes(x = h_vector, y = value, color = column)) +
plot_lines +
labs(title = "Circle graph",
x = TeX("$\\hat{h}$"),
y = "Covariance Error",
color = expression(beta)) +
scale_x_log10(breaks = h_vector, labels = TeX(h_label_latex)) +
scale_y_log10() +
theme_minimal() +
theme(text = element_text(family = "Palatino"))